Experimental  Implementation  of  a  Model-Based 
Inverse  Filter  to  Attenuate  Hysteresis  in  an 
Atomic  Force  Microscope 

Andrew  Hatch*,  Ralph  C.  Smith*  and  Tathagata  De t 

*Center  for  Research  in  Scientific  Computation,  North  Carolina  State  University,  Raleigh,  NC  27695 
^Electrical  Engineering  Department,  Iowa  State  University,  Ames,  IA  50011 


Abstract 

This  paper  addresses  the  development  and  exper¬ 
imental  validation  of  a  model-based,  open  loop  con¬ 
trol  design  for  mitigating  the  frequency-dependent  ef¬ 
fects  of  hysteresis  in  an  atomic  force  microscope  (AFM). 
The  models  are  based  on  homogenized  energy  relations 
which  characterize  the  hysteretic  constitutive  behavior 
of  the  piezoceramic  AFM  stage.  Approximate  model 
inverses  are  then  employed  as  filters  to  linearize  trans¬ 
ducer  dynamics  for  control  design.  When  experimen¬ 
tally  implemented  in  open  loop  control  designs,  inverse 
compensation  in  this  manner  produces  an  approximately 
tenfold  increase  in  tracking  accuracy  as  compared  with 
the  unfiltered  case. 


1.  Introduction 

Crucial  components  in  atomic  force  microscope  (AFM) 
designs  are  the  piezoceramic  (PZT)  drive  mechanisms 
used  to  position  the  sample  [3,  9].  Whereas  the  spe¬ 
cific  designs  vary,  all  employ  the  converse  piezoelectric 
effect  to  provide  input  displacements  and  force  in  re¬ 
sponse  to  applied  voltages  or  fields.  This  provides  mech¬ 
anisms  which  have  nanometer-level  setpoint  accuracy 
but  comes  at  the  cost  of  the  hysteresis  and  constitutive 
nonlinearities  inherent  to  piezoelectric  compounds. 

To  illustrate,  consider  the  prototypical  stage  depicted 
in  Figure  1(a)  which  employs  stacked  piezoceramic  ac¬ 
tuators  to  position  the  sample  in  the  x  and  y  directions. 
An  additional  mechanism  provides  transverse  position¬ 
ing  capabilities.  Nested  minor  loops  collected  at  0.1 
Hz  are  plotted  in  Figure  1(b)  and  data  collected  at  fre¬ 
quencies  ranging  from  0.279  Hz  to  27.9  Hz  is  plotted  in 
Figure  2  to  illustrate  the  frequency-dependent  nature  of 
the  hysteresis. 

At  low  frequencies,  the  inherent  hysteresis  can  be 
accommodated  through  PID  or  robust  control  designs 
[1,  2,  4,  5,  7].  However,  at  the  higher  frequencies  re¬ 
quired  for  applications  ranging  from  real-time  monitor¬ 
ing  of  biological  processes  —  e.g.,  protein  unfolding 
to  comprehensive  product  diagnostics,  increasing  noise- 
to-data  ratios  and  diminishing  high-pass  characteristics 
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Figure  1.  (a)  PZT-based  AFM  stage,  (b)  Nested  minor 
loops  in  AFM  data  collected  at  0.1  Hz. 

of  control  filters  preclude  a  sole  reliance  on  feedback  laws 
to  eliminate  hysteresis.  This  motivates  the  development 
of  control  designs  that  incorporate  and  approximately 
compensate  for  hysteresis  through  model  inverses  em¬ 
ployed  either  in  feedback  or  feedforward  loops. 

In  this  paper,  we  employ  previously  developed  ho¬ 
mogenized  energy  models  to  characterize  the  hysteretic 
relation  between  input  fields  and  strains  generated  by 
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Figure  2.  Frequency-dependent  field-displacement  behavior  of  an  AFM.  Sample  rates  of  (a)  0.279  Hz,  (b)  5.58  Hz, 
and  (c)  27.9  Hz. 


the  PZT  rods.  These  constitutive  laws  are  then  em¬ 
ployed  to  construct  inverse  models  which  serve  to  lin¬ 
earize  the  transducer  response.  When  experimentally 
implemented  in  open  loop  control  designs,  it  is  illus¬ 
trated  that  inverse  compensation  in  this  manner  yields 
an  approximately  tenfold  increase  in  accuracy  over  the 
unfiltered  case  when  operating  at  27.9  Hz. 


minimum  of  t/j  occurs.  The  parameter  r/  is  the  reciprocal 
of  the  slope  of  the  E-P  relation  after  switching  occurs. 
This  fact  can  be  used  to  establish  an  initial  parameter 
value  for  when  modeling  a  specific  data  set. 

In  the  case  of  no  applied  stresses  a,  the  Gibbs  energy 
can  be  formulated  as 

G(E,  P)  =  ip  —  EP  (2) 


2.  Constitutive  Relations 

To  model  the  constitutive  behavior  of  the  piezoce¬ 
ramic  stacked  actuator,  the  stress-strain  relation  is  as¬ 
sumed  to  be  linear.  However,  the  relation  between  the 
applied  voltage  V  (or  the  applied  field  E)  and  the  po¬ 
larization  P  exhibits  nonlinearities  and  hysteresis.  The 
actuator  is  also  assumed  to  be  biased  through  poling  so 
that  the  relation  between  P  and  the  strain  e  is  linear 
for  the  considered  operating  conditions.  To  character¬ 
ize  the  hysteretic  E-P  behavior  at  the  domain  level, 
a  Helmholtz  energy  relation  was  derived  in  [10]  using 
statistical  mechanics  principles  under  the  assumption 
that  dipoles  are  either  aligned  with  the  field  or  diamet¬ 
rically  opposed  to  it.  This  model  is  appropriate  for  a 
single  crystal  with  uniform  effective  fields.  To  construct 
a  macroscopic  model  for  a  polycrystalline  material  with 
variable  effective  fields,  the  coercive  and  effective  field 
values  are  then  assumed  to  be  distributed  as  detailed  in 
the  latter  half  of  this  section. 

Under  fixed  temperature  conditions  with  no  applied 
stress  tr,  it  is  illustrated  in  [10]  that  a  first  order  approxi¬ 
mation  to  the  statistical  mechanics-based  Helmholtz  en¬ 
ergy  is  the  piecewise  quadratic  relation 

P  +  Pr)2  P<-Pi 

m  =  \  Hp  -  pn)2  p^pi  (1) 

[iriiP!  -  PR)  (g  -  PR)  \P\<Pi. 

As  shown  in  Figure  3,  Pi  is  the  positive  inflection  point 
and  Pr  is  the  value  of  P  at  which  the  positive  local 


where  the  second  term  represents  the  electrostatic  en¬ 
ergy  due  the  applied  field  E.  In  order  to  include  fer- 
roelastic  coupling,  we  utilize  the  extended  Helmholtz 
relation 


MP,  e)  =  V>( P )  +  \yPe2  -  UP7 eP.  (3) 

The  Gibbs  energy  is  then  given  by 
G(E,  P,  e)  =  V>(P)  +  *FPe2  -  YP'yeP  —  EP  —  ae  (4) 


Figure  3.  (a)  Helmholtz  energy  iji  and  Gibbs  energy  G 
for  (7  =  0  and  increasing  fields  E.  (b)  Dependence  of  the 
local  polarization  P  on  the  field  E  at  the  lattice  level  in 
the  absence  of  thermal  activation. 
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where  ae  incorporates  the  elastic  energy.  Note  that  Yp 
is  the  Young’s  modulus  for  a  constant  polarization  and 
7  is  a  ferroelastic  coupling  coefficient. 

In  the  case  of  negligible  thermal  activation,  the  local 
average  polarization  P  is  determined  from  the  necessary 
conditions 


dG 

dP 


=  0 


d2G 

dP 2 


>  0. 


(5) 


Applying  these  conditions  to  (2)  yields  a  piecewise  linear 
E-P  characterization 


[  [P(£,£;£c,£)](0) 

[P(E,e;Ec,Zm=  \  f  "  (6) 

E  ,  PrV 
l.  V  '  rj—2Ypje 

for  the  respective  cases  (r(i)  =  0},  {r(t)  ^  0  and 
P(maxr(t))  =  —  Ec},  {r(t)  ^  0  and  P(maxr(t))  = 
Ec}.  The  transition  points  r  are  specified  by 


r{t)  =  {t  G  (0,  Tf]  |  E(t)  =  -Ec  or  E(t)  =  Ec} , 

and 

{E  _  Pay 
r i  r\—2Yp^e 

i 

E  ,  PrV 

V  '  /?  — 21'p7t 

denotes  the  initial  dipole  orientation  for  respective  ini¬ 
tial  fields  {P(0)  <  —  Ec},  {—Ec  <  P(0)  <  Ec}  or 
{P(0)  >  Ec}. 

However,  if  thermal  activation  is  significant,  dipoles 
can  achieve  the  thermal  energy  required  to  switch  in  ad¬ 
vance  of  the  minimum  Gibbs  energy  so  the  relative  ther¬ 
mal  and  Gibbs  energy  must  be  balanced  through  Boltz¬ 
mann  principles.  The  probability  density  for  achieving 
an  energy  level  G  is  then  given  by 

H(G)  =  Ce~GV/kT  (7) 


where  k  is  Boltzmann’s  constant,  V  is  a  reference  vol¬ 
ume  and  C  is  a  constant  that  is  selected  so  that  when 
I-l(G)  is  integrated  over  all  possible  dipole  orientations, 
a  probability  of  1  is  achieved.  If  we  let  2cr  be  the  sep¬ 
aration  between  possible  polarization  states  around  Po, 
the  probabilities  of  reaching  a  polarization  state  having 
sufficient  energy  to  switch  orientations  are  given  by 


r+-  = 


('PO+CT 

>P0-<?  1 


/•OO 
Pn-<T  ' 


r_+  = 


rPo+cr  _ 

JPo-a  e 


f- 


Po+& 


-G(E,P)V/kTdp 

-G(E,P)V/kTdp 

G(E,P)V/kTdp 
G(E,P)V/kTdp  ' 


(8) 

(9) 


The  likelihoods  of  reaching  the  required  energy  and  thus 
of  the  dipoles  switching  from  a  positive  to  a  negative 
orientation  and  conversely  are  then 

P+-  =  ~r+-  ,  P-+  =  ~r-+  (10) 

r  t 


where  r  is  the  relaxation  time.  The  fractions  of  dipoles 
in  each  orientation  evolve  according  to  the  ordinary  dif¬ 
ferential  equations 

dx+ 

—  =  -P+-X++P-  +  X- 
^  =  -P-+x-  +P+-X+. 


The  expected  polarizations  due  to  positively  and  neg¬ 
atively  oriented  dipoles  are 

<•00  rPo-a 

(P+>  =  /  Pp(G)dP  ,  <P+)  =  /  Pp{G)dP 

J  Po+CT  J  —OO 

so  that  evaluation  of  C  yields 

f~  Pe-G(E,P,T)V/kTdp 
/p  \  JPo+<7 

'  +/  f00  g -G(E,P,T)V/kTdP 

JP0+<T 

(11) 

/•Po-o-  Pe-G{E,P,T)V/kTdp 

(P  )  =  J~°° 

[p0-<?  e-G(E,P,T)V/kTdP 

(12) 

For  a  single  crystal  with  uniform  effective  field, 
average  polarization  is  subsequently 

the  local 

P  =  X+{P+)  +  X-(P-). 

(13) 

In  the  manner  detailed  in  [10],  the  evaluation  of  the 
integrals  in  (8)  and  (11)  can  be  simplified  through  ap¬ 
proximations  employing  the  inflection  points  ±P/  rather 
than  the  unstable  equilibrium  Pq. 

Both  of  the  relations  (6)  and  (13)  are  valid  only  for 
homogeneous  single  crystal  materials  with  uniform  ef¬ 
fective  fields.  To  account  for  nonuniformity  and  inho¬ 
mogeneities  in  the  materials,  local  coercive  and  effective 
fields  are  assumed  to  be  manifestations  of  underlying 
distributions  rather  than  constants.  The  macroscopic 
polarization  model  is  then  given  by 

POO  pOO 

P{E)=  /  /  P(E  +  Ee,Ec,()u1{Ec)u2(Ee)dEedEc 

J  0  J  —  oo 

(M) 

where  iq  and  v2  are  appropriate  densities.  Motivated  by 
choices  in  the  magnetics  literature,  zq  and  v2  were  re¬ 
spectively  designated  to  be  lognormal  and  normal  den¬ 
sities  in  [10].  However,  the  fact  that  neither  of  these 
choices  is  based  on  energy  considerations,  motivates  con¬ 
sideration  of  general  densities. 

Based  on  physical  arguments,  we  may  assume  that 
the  general  densities  decay  to  zero.  Therefore,  the  dou¬ 
ble  integral  can  be  evaluated  by  truncating  the  domains 
and  using  composite  Gaussian  quadrature.  In  this  case, 
the  discretization  of  (14)  gives 

Ni  Nj 

P{E)  =  y^y  P{E  +  Ee.,ECi,£i)v\{Ee.)v2(Ec.)viWj 

i= 1  7  —  1 

(15) 
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where  Eej ,  ECi  are  the  abscissas  and  t>, ,  Wj  are  the  weights 
Here,  v\(Ee.)  and  ^2 ( ECt )  are  parameter  values  to  be  de¬ 
termined.  Details  regarding  the  identification  of  v\  and 
V2  can  be  found  in  [6]. 


Finally,  the  equilibrium  condition 

os 

(16) 

yields  the  elastic  constitutive  relation 

a  =  Ype  -  YpjP. 

(17) 

where  the  resultant  J\f  =  fA  crdA  is  given  by 

Af  =  YpAp  +  cpA-f^--YpA7P(E).  (19) 
ox  ox  at 

Note  that  u  denotes  displacement  in  the  longitudinal  re¬ 
direction  and  that  the  relation  strain  £  =  §7  is  used  in 
obtaining  (19).  The  boundary  conditions  are  u(t,  0)  =  0 
at  the  fixed  end  and 

Ou  d2u 

A f(tj)  =  -kLu(t,£)  -  -  ML^(t,£)  (20) 


This  relation  along  with  the  nonlinear  polarization  re¬ 
lation  (14)  quantifies  the  constitutive  behavior  for  the 
piezoceramic  materials  employed  in  the  AFM. 


3.  Model  Inverse 

To  construct  a  model  inverse  to  be  employed  in  sub¬ 
sequent  control  design,  we  exploit  the  approximate  mono¬ 
tonicity  of  the  E-P  relation  to  construct  an  algorithm 
in  which  the  forward  model  is  advanced  until  the  pre¬ 
scribed  polarization  is  crossed.  Interpolation  is  then 
used  to  specify  a  final  polarization  value  corresponding 
to  the  input  field.  This  process  is  outlined  Algorithm  1. 

Algorithm  1. 

for  k  =  2  :  Nk 
Specify  A E 

dP  =  Pk~  Pk- 1 

Etmp  —  Ek  —  1  5  Ptmp  =  Pk—1 

while  sgn (dP)  ■  (Pk  -  Pk- 1)  >=  0 
Etmp  =  Etmp  T  A E 
Ptmp  given  by  (15) 

end 

Ek  given  by  linear  interpolation 

end 


4.  Stacked  Actuator  Model 


In  addition  to  hysteresis,  the  dynamics  of  the  actua¬ 
tor  must  be  incorporated.  We  assume  that  the  rod  has 
cross-sectional  area  A,  length  £,  density  p  and  Young’s 
modulus  Yp .  Let  cp  be  the  Kelvin-Voigt  damping  pa¬ 
rameter  and  7  be  the  piezoelectric  coupling  coefficient. 
In  the  present  stage  design,  depicted  in  Figure  1(a),  one 
end  of  the  actuator  is  fixed,  while  the  attachment  at 
the  other  end  can  be  modeled  as  a  damped  spring-mass 
system.  For  this  end,  let  ML  be  the  mass,  be  the 
stiffness,  and  cl  be  the  damping  coefficient.  Force  bal¬ 
ancing  along  the  actuator  then  yields  the  relation 


d2u  _  (W 

pAW  ~  foT 


(18) 


at  the  moving  end.  Initial  conditions  are  given  by  u( 0,  x)  = 
=  0.  The  polarization  P{E)  is  specified  by  (14). 

An  approximate  solution  to  (18)  is  found  by  first 
deriving  a  weak  form 


fe  „  d2u 

fe 

0  ,  du 

p  92u 

i  fAaWHx  +  , 

L 

Y PA— 

OX 

+  C  Adxdt 

-r—dx 

ox 


=  -  [  YpA(3P(E)^-dx 
Jo  dx 

(2!) 

and  then  discretizing  in  space  with  linear  finite  elements 
and  in  time  with  finite  difference  techniques.  Details  of 
the  derivation  of  the  weak  form  and  the  construction  of 
the  finite  element  and  finite  difference  equations  can  be 
found  in  [9]. 

In  order  to  determine  the  displacement  of  the  rod 
at  the  end  x  =  £,  the  rod  model  can  be  approximated 
by  a  lumped  spring-mass  model  since  the  cross-section 
is  small  compared  to  the  length  and  the  electric  field  E 
is  approximately  uniform  throughout  the  actuator.  Nu¬ 
merical  simulations  and  physical  experiments  demon¬ 
strate  that  the  resulting  ODE  model  is  a  good  approx¬ 
imation  to  the  PDE  model  for  specifying  tip  displace¬ 
ment.  The  ODE  model  has  the  advantage  of  faster  com¬ 
putational  times  compared  to  the  discretized  PDE. 


5.  Model  Validation 

To  demonstrate  the  accuracy  and  efficiency  of  the 
hysteresis  model,  we  consider  the  characterization  of 
nested  minor  loops  collected  at  0.1  Hz  as  well  as  E-P  be¬ 
havior  at  frequencies  ranging  from  0.279  Hz  to  27.9  Hz. 
In  both  cases,  displacements  were  computed  using  the 
ODE  model  constructed  using  the  stress  relation  (17) 
and  hysteresis  model  (14)  with  general  densities  v\  and 
v-2  identified  using  the  techniques  detailed  in  [6,  8]. 

Figure  4  demonstrates  the  characterization  of  fre¬ 
quency  dependent  dynamics.  The  latter  involves  the 
quantification  of  both  thermal  relaxation  and  inertial 
effects  as  illustrated  by  the  change  in  sign  of  the  slope 
QE  following  field  reversal.  Further  details  demonstrat¬ 
ing  properties  of  the  model  for  characterizing  hysteresis 
in  various  PZT  compounds  can  be  found  in  [6,  8,  10]. 


4 


Figure  4.  Characterization  of  AFM  field-displacement  behavior  with  sample  rates  of  (a)  0.279  Hz,  (b)  5.58  Hz 
and  (c)  27.9  Hz. 


6.  Open  Loop  Control  Implementation 

To  illustrate  the  effect  of  filters  employing  the  in¬ 
verse  model  from  Section  3  on  open  loop  tracking  perfor¬ 
mance,  we  summarize  experiments  conducted  at  0.279  Hz 
and  27.9  Hz.  In  both  cases,  the  specified  trajectory  to 
be  tracked  consisted  of  a  triangle  wave  ranging  from  0  to 
7000  V/m.  To  construct  the  model-based  inverse  filter, 
the  specified  trajectory  was  input  to  the  inverse  model 
and  the  resulting  field  was  applied  to  the  experimental 
device.  To  provide  a  metric  for  comparison,  a  second 
predicted  input  field  for  each  case  was  determined  by  a 
linear  scaling  of  the  desired  output  displacement.  The 
scaling  factors  were  derived  from  manufacturer  specifi¬ 
cations  for  the  actuator  stroke  when  the  field  is  increased 
throughout  its  range. 

The  specified  and  achieved  trajectories,  errors,  and 
hysteresis  plots  for  the  two  frequencies  are  plotted  in 
Figures  5  and  6.  At  0.279  Hz,  the  filtered  design  pro¬ 
vides  only  marginally  improved  accuracy  due  to  the  low 
degree  of  hysteresis  and  constitutive  nonlinearities.  For 
the  more  hysteretic  response  at  27.9  Hz,  however,  the 
inverse  filter  provides  a  significant  increase  in  accuracy 
and  yields  errors  that  are  approximately  a  factor  of  10 
smaller  than  the  unfiltered  but  scaled  case.  This  il¬ 
lustrates  the  advantage  of  incorporating  the  frequency- 
dependent  model  inverse  in  the  control  design. 

We  note  that  the  primary  source  of  errors  in  the 
filtered  design  is  variability  between  experiments  as  il¬ 
lustrated  by  the  variation  in  the  hysteresis  plots  at  each 
frequency.  This  is  hypothesized  to  be  due  to  variations 
in  the  true  applied  voltage  and  illustrates  one  reason 
feedback  is  necessary  in  final  control  designs. 


7.  Concluding  Remarks 
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The  modeling  framework  employed  in  this  paper  pro-  Figure  5.  Tracking  performance  utilizing  the  model- 

vides  the  capability  for  characterizing  the  frequency-  based  filter  and  without  a  filter  at  0.279  Hz. 

dependent  hysteresis  inherent  to  the  piezoceramic  drive 
mechanisms  in  atomic  force  microscopes.  It  is  also  ame- 
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